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Abstract 



We determine the critical equation of state of three-dimensional randomly 
dilute Ising systems, i.e. of the random-exchange Ising universality class. We 
first consider the small-magnetization expansion of the Helmholtz free energy 
in the high-temperature phase. Then, we apply a systematic approximation 
scheme of the equation of state in the whole critical regime, that is based on 
polynomial parametric representations matching the small-magnetization of 
the Helmholtz free energy and satisfying a global stationarity condition. These 
results allow us to estimate several universal amplitude ratios, such as the ratio 
A + /A~ of the specific-heat amplitudes. Our best estimate A + /A~ = 1.6(3) 
is in good agreement with experimental results on dilute uniaxial antiferro- 
magnets. 



PACS Numbers: 75.10.Nr, 64.60.Ak, 75.10.Hk 
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I. INTRODUCTION. 



The critical properties of randomly dilute Ising systems have been much investigated 
experimentally and theoretically, see, e.g., Refs. [1-5] for reviews. The typical example is 
the ferromagnetic random Ising model (RIM) with Hamiltonian 

7Y= -J PiPj SiSj - HY^PiSi, (1.1) 

<ij> i 

where J > 0, the first sum extends over all nearest-neighbor sites, S; = ±1 are Ising 
spin variables, and pi are uncorrelated quenched random variables, which are equal to one 
with probability p (the spin concentration) and zero with probability 1 — p (the impurity 
concentration). Considerable work has been dedicated to the identification of the universal 
critical behavior of the RIM, and in particular to the determination of the critical exponents, 
which have been computed with great accuracy in experiments and in theoretical works. On 
the other hand, much less is known about the critical equation of state and the corresponding 
universal amplitude ratios. From the experimental side, this is essentially due to the fact 
that typical experimental realizations of the RIM are uniaxial antiferromagnets such as 
Fe x Zn 1 _ ;c F2 and Mn x Zn 1 _ ;I .F2 materials, which are usually modeled by the Hamiltonian 
(1.1) with J < 0. For H = 0, there is a simple mapping between the ferromagnetic and 
the antiferromagnetic model, so that the critical behavior of the antiferromagnets can still 
be obtained from that of the ferromagnetic RIM. The situation is more complex for H ^ 0. 
The RIM (1.1) with H ^ corresponds to an antiferromagnet with a staggered magnetic 
field that cannot be realized experimentally. Conversely, antiferromagnets in a uniform 
magnetic field have a different critical behavior and belong to the same universality class of 
the ferromagnetic random-field Ising model [6]. For these reasons the equation of state of 
the RIM is not relevant for dilute uniaxial antiferromagnets. However, from the equation 
of state one can derive amplitude ratios involving quantities defined in the high- and low- 
temperature phase that can be measured experimentally, see, e.g., Ref. [3]. For example, 
the ratio A + /A~ of the specific- heat amplitudes in the high- and low-temperature phase has 
been determined: A+/A~ = 1.6(3) (Ref. [7]), and A+/A~ = 1.55(15) (Ref. [3]). On the 
theoretical side, only a few works, based on field-theoretical (FT) perturbative expansions, 
have attempted to determine other universal quantities beside the critical exponents [8-12], 
obtaining rather imprecise results. For example, these calculation have not even been able 
to determine reliably the sign of the ratio A + /A~ . Actually, e-expansion calculations [9,10] 
favor a negative value, in clear disagreement with experiments. 

In this paper we determine the equation of state of the RIM in the critical region, that 
is the relation among the external field H, the reduced temperature t = (T — T c )/T c , and 
the magnetization 

M = ^(ZPi*i)> (1-2) 

i 

where the overline indicates the average over the random variables pi, and ( ) indicates the 
sample average at fixed disorder. It can be written in the usual form as 

H oc M s f(x), x oc tM 1//3 , (1.3) 
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where x and f(x) are normalized so that x — — 1 corresponds to the coexistence curve, hence 
/(— 1) = 0, and /(0) = 1. In order to determine the critical equation of state in the whole 
critical region, we use an approximation scheme that has already been applied with success 
to the pure Ising model in three [13-15] and in two dimensions [16]. We first consider the 
expansion of the equation of state in terms of the magnetization in the high-temperature 
phase. The first few nontrivial coefficients can be determined either from Monte Carlo 
simulations of the RIM or from the analysis of FT perturbative expansions. These results 
are then used to construct approximations that are valid in the whole critical region and 
that allow us to determine several universal amplitude ratios. For example, we anticipate 
our best estimate of the specific-heat amplitude ratio 



which compares very well with the experimental determinations in dilute uniaxial antiferro- 
magnets [7,3]. 

The paper is organized as follows. In Sec. II we discuss the general properties of the 
critical equation of state. In Sec. Ill we consider its small-magnetization expansion in the 
high-temperature phase, reporting estimates of the first few nontrivial terms, obtained by 
Monte Carlo simulations and FT methods. These results are used in Sec. IV to construct 
approximate polynomial parametric representations that are valid in the whole critical region 
and that allow us to achieve a rather accurate determination of the scaling function f(x). 
Finally, in Sec. V we determine several universal amplitude ratios, such as the specific-heat 
amplitude ratio A + / A~ . In App. A we report the definitions of the thermodynamic quantities 
that are considered in the paper. In App. B we discuss the correspondence between the RIM 
correlation functions and the correlation functions of the corresponding translation-invariant 
field theory that is obtained by using the standard replica trick. App. C reports some details 
on the six- loop FT calculation of the universal amplitude ratio Rt . 



The equation of state relates the magnetization M, the magnetic field H, and the reduced 
temperature t = (T — T c )/T c . In the neighborhood of the critical point t — 0, H — 0, it can 
be written in the scaling form 



where B c and B are the amplitudes of the magnetization on the critical isotherm and on 
the coexistence curve, see App. A. According to these normalizations, the coexistence curve 
corresponds to x — —1, and the universal function f(x) satisfies /(— 1) = and /(0) = 1. 
The apparently most precise estimates of the critical exponents have been recently obtained 
by Monte Carlo simulations of the RIM: Ref. [17] obtains v = 0.683(3) and rj = 0.035(2), 
while Ref. [18] reports v = 0.6837(53) and rj = 0.0374(45). The FT estimates v = 0.678(10) 
and rj = 0.030(3), obtained by analyzing six-loop fixed-dimension series [19], are in good 



A+/A- = 1.6(3), 



(1.4) 



II. THE CRITICAL EQUATION OF STATE 



H = B; s M s f(x), 
x = t(M/B)- 1/l3 , 



(2.1) 
(2.2) 
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agreement [20]. The other exponents a, 7, f3, and 5 can be determined using scaling and 
hyperscaling relations. 

The equation of state is analytic for \H\ > 0, implying that f(x) is regular everywhere 
for x > — 1. In particular, f(x) has a regular expansion in powers of x around x = 

00 

/(*) = 1 + E (2-3) 

n=l 

At the coexistence curve, i.e., for x — > — 1, /(x) is expected to have an essential singularity 
[21], so that it can be asymptotically expanded as 

00 

/(^E/r x (i+*) n (2.4) 

n=l 

The free energy of a dilute ferromagnet is nonanalytic at H = for all temperatures below 
the transition temperature of the pure system [22] , which is larger than the critical temper- 
ature of the dilute one. Therefore, at variance with pure systems, the equation of state of 
the RIM is also nonanalytic for t > and H = 0. However, as argued in Ref. [23], these 
singularities are very weak and all derivatives remain finite at H = 0, as it happens at the 
coexistence curve. This allows us to write down an asymptotic large- x expansion of the form 

00 

/(*) = x-r £ f^x-^. (2.5) 

n=0 

It is useful to rewrite the equation of state in terms of a variable proportional to Mt^ 13 , 
although in this case we must distinguish between t > and t < 0. For t > we define 

1/2 



/r+\ L/z 



ct 



(c+y 



n l/2 

Mr", (2.6) 



while for £ < we set 



- - f H)-'- (2-7) 

The constants C 1 * 1 and are the amplitudes appearing in the critical behavior of the 

two- and four-point susceptibilities x an d X4 respectively, see App. A. With the chosen 
normalizations [4], 

F W= , + I,3 + g_l_ r24 ^- 1 , (2 . 8) 

OO 1 

*(«) = (« - 1) + E 7— w^- ( M - ( 2 - 9 ) 



J=3 
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The large- z expansion of the scaling function F(z) is given by 

F(z) = z s J2 F™z- k/f3 . (2.10) 

fc=0 

The functions F(z) and $(w) are clearly related to f(x). Indeed, 

z~ 5 F(z) = F °°f(x), z = z x-P, (2.11) 

«"'*(«) = ^^-m, u = (-x)-P. (2.12) 

where z = (Rt) 1 ^ 2 is a universal constant, see Sec. V. 

In order to determine the critical equation of state, we first consider its small- 
magnetization expansion. Then, we construct parametric representations of the critical 
equation of state based on polynomial approximations, which are valid in the whole criti- 
cal region. This method have already been applied to the Ising universality class in three 
[13-15] and two dimensions [16], and to the three-dimensional XY and Heisenberg univer- 
sality classes [24,25]. 

III. SMALL-MAGNETIZATION EXPANSION IN THE HIGH-TEMPERATURE 

PHASE 

A. Small- magnetization expansion of the Helmholtz free energy 

In the high-temperature phase the quenched Helmholtz free energy A(t, M) admits an 
expansion around M = 0: 

e [A(t, M) - Ait, 0)] = l -M 2 + E J^yG 23 M 2 \ (3.1) 

where £ is the second-moment correlation length along the axis // M 0. M r.U£ 
c = (/ + ) 1-r?/ ' 2 /(C + ) 1 / 2 , and C + , / + are the amplitudes of x an d £ respectively, see App. A 
for notations. We recall that the equation of state is related to the Helmholtz free energy 
by 

dA(t, M) 

and therefore the expansion (3.1) is strictly related to the expansion of the function F(z), 
cf. Eq. (2.6). 

In the critical limit the coefficients G 2 j of the expansion (3.1) are universal. They can 
be determined from the high-temperature critical limit of combinations of zero-momentum 
connected correlation functions averaged over the random dilution 



d n ~ l M 

Xn 



dH 



n-l 



(3.3) 

H=0 
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Indeed, 



G 4 = lim lim (--^) , (3.4) 
G G = lim lim I --^ + Xa ' 



G 8 = lim lim + - 280- A { 



where %2 = X and V is the volume. 

For the purpose of determining the small-magnetization expansion of the RIM equation 
of state, it is convenient to rewrite the Helmholtz free energy of the RIM, cf. Eq. (3.1), in 
the equivalent form 

A(t, M) - A(t, 0) = J^t^A(z), (3.5) 
where z is defined in Eq. (2.6), 



and 

- / > (3-7) 



Since F(z) = dA(z)/dz, Eq. (2.6) follows from Eq. (3.6). 

Some of these universal constants have been recently estimated by means of a Monte 
Carlo simulation [17], obtaining 

G 4 = 43.3(2), r 6 = 0.90(15). (3.8) 

Correspondingly, G 6 = 1.7(3) x 10 3 . 



B. Field-theoretical approach 

One may estimate the universal quantities Gij and r2j by FT methods. The FT approach 
is based on an effective Landau-Ginzburg- Wilson Hamiltonian [26] that can be obtained by 
using the replica method, 

r i 1 N 1 N N 1 

?V = J d 3 x - ]T [(d^ t ) 2 + ry^ 2 ] +-JT («o + ) + ^ & ' ( 3 - 9 ) 

where ipi is an iV-component field. The critical behavior of the RIM is expected to be 
described by the Hamiltonian for u < and in the limit iV — > 0. Conventional FT 
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computational schemes show that the fixed point corresponding to the pure Ising model is 
unstable and that the renormalization-group flow moves towards another stable fixed point 
describing the critical behavior of the RIM. 

The most precise FT results for the critical exponents have been obtained in the frame- 
work of the perturbative fixed-dimension expansion in terms of zero-momentum couplings. 
The corresponding Callan-Symanzik /3-functions and the renormalization-group functions 
associated with the critical exponents have been computed to six loops [19,27]. 

The Helmholtz free energy A^ associated with the Hamiltonian (3.9) can be written as 
[28,29] 

e iA^(M) - A^(0)] = \H^l + T\ 5>4i + 9A25 ab )M 2 a M 2 b + (3.10) 

a ' ab 

+4 E(<?6i + 9^ab + g e3 5 ahc )M 2 M 2 M 2 + 

abc 

+ ^ E (081 + 982$ab + gssSabScd + <? 8 Abc + <?85<Ucd) M'M^Mj + ... 
abed 

where M a = cM a ^/ u , c = ^ 1 ~ 11 ^ 2 /x 1 ^ 2 \h=o- Note that g 4 \ = u* and g 42 = v*, where u* and 
v* are the fixed-point values of the renormalized quartic couplings u and v. Applying the 
replica method in the presence of a uniform external field H, see App. B for details, one can 
identify 

G 4 = g 42 = v*, G 6 =g 63 , G s = g S 5, (3.11) 

etc... It is worth noting that the other coefficients appearing in the small- magnetization 
expansion of the Helmholtz free energy (3.10), i.e. g 4 i, gei, g$2, etc., can be related to 
dilution averages of products of single-sample n-point correlation functions, as shown in 
App. B. 

The quartic coupling G 4 can be estimated from the position of the RIM fixed point in 
the u, v plane, i.e. from the common zero of their /5-functions. The analysis of the six-loop 
perturbative series gives results somewhat dependent on the resummation method [19], see 
also the five-loop analysis of Ref. [30]. Combining all results together, one finds 36.5 ^$ G 4 ^ 
39.5, which can be summarized as G 4 = 38.0(1.5). Such a result is not consistent with the 
Monte Carlo estimate (3.8), obtained by estimating the limit given in Eq. (3.4). But, as 
discussed in Ref. [17], the quantitative consequences for the determination of the critical 
exponents are very small, since the renormalization-group functions corresponding to the 
critical exponents turn out to be little sensitive to the correct position of the fixed point [31] 
(along the Ising-to-RIM renormalization-group trajectory, see also Ref. [32]). 

The universal constants G 6 and r 6 = G 6 /G 4 have been estimated in Ref. [29] by analyzing 
the corresponding four-loop series with the Pade-Borel method. These analyses provide the 
estimate r 6 fa 1.09. We have reanalyzed this series, finding that, unlike critical exponents, 
tq is rather sensitive to the position of the fixed point. Using the Monte Carlo estimates of 
u* and v* and the Pade-Borel method, we obtain r$ fa 0.6. Thus, a conservative estimate 
is r-Q = I.IIq^, which is in agreement with the Monte Carlo result of Ref. [17] mentioned 
in Sec. Ill A, r 6 = 0.90(15). We have also computed the fixed-dimension perturbative 
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expansion of r 8 = G&/G\ to three loops. This requires the computation of the 8-point one- 
particle irreducible zero-momentum correlation function. At three loops this calculation 
requires the evaluation of 42 Feynman diagrams. For this purpose we have used the general 
algorithm of Ref . [33] . The expansion of r 8 is 

105 35 
r 8 = - — (8u + 3v) + ——(2480m 2 + 2286wv + 585-u 2 ) (3.12) 
6471 92167T 2 

-0.132874m 3 - 0.199366 u 2 v - 0.1080052 uv 2 - 0.0211475 v 3 + ... 

where the zero-momentum quartic couplings u and v are normalized so that u = uq/itl 
and v = Vo/m at tree level (see App. C for precise definitions). We have checked that 
Eq. (3.12) correctly reproduces the known series for the O(N) models [34] with N — 0,1 
in the appropriate limits. Unfortunately, the analysis of the expansion (3.12) provides only 
a very rough estimate [35], r 8 = 15(15). We will obtain a much better estimate of r 8 in 
Sec. IV B, from the results for the equation of state. 



IV. APPROXIMATE REPRESENTATIONS OF THE CRITICAL EQUATION OF 

STATE 

A. Polynomial parametric representations 

In order to obtain approximate expressions for the equation of state, we parametrize the 
thermodynamic variables in terms of two parameters R and 9, implementing all expected 
scaling and analytic properties. Explicitly, we write [36] 

M = moR^e, 

t = R{i-e 2 ), 

H = h R /3S h{9), (4.1) 

where ho and m are normalization constants. The variable R is nonnegative and measures 
the distance from the critical point in the (t, H) plane; the critical behavior is obtained for 
R — > 0. The variable parametrizes the displacements along the lines of constant R. The 
line 9 = corresponds to the high-temperature phase t > and H = 0; the line 9 = 1 to 
the critical isotherm t = 0; 9 = 9 , where 9q is the smallest positive zero of h{9), to the 
coexistence curve T < T c and H — > 0. Of course, one should have 9 > 1 and h(9) > 
for < 9 < 9q. The function h(9) must be analytic in the physical interval < 9 < 9 in 
order to satisfy the requirements of regularity of the equation of state (Griffiths' analyticity). 
Note that the mapping (4.1) is not invertible when its Jacobian vanishes [4], which occurs 
for 9f = 1/(1 — 2(3). Thus, a parametric representation is acceptable only if 9 < 9\. The 
function h{9) must be odd in 9, to guarantee that the equation of state has an expansion 
in odd powers of \M\ in the high-temperature phase for \M\ — > 0. Moreover, it can be 
normalized so that h{9) = 9 + 0(9 3 ). 

The scaling functions f(x) and F(z) can be expressed in terms of 9. The scaling function 
f(x) is obtained from 
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X 



e 2 (e»\ 1/p 



- S h{6) 



while F(z) is obtained by 



z = pe(i-e*y\ 

F(z(e)) = P (i-e*y p5 h(e), (4.3) 

where p can be related to m , ho, C + , and CX by using Eqs. (2.6) and (4.1). 

Eq. (4.1) and the normalization condition h(9) m 9 for 9 — > do not completely fix the 
function h(9). Indeed, one can rewrite the relation between x and 9 in the form 

x y = h(l) / °°(1 - 9 2 y9 1 - 5 . (4.4) 

Thus, given f(x), the value of h(l) can be arbitrarily chosen to completely fix h(9). One 
may fix this arbitrariness by choosing arbitrarily the parameter p in the expression (4.3). 
We approximate h{9) with polynomials, i.e., we set 

h(9)=9+j2h 2n+1 9 2n+1 . (4.5) 

n=l 

This approximation scheme turned out to be effective in the case of pure Ising systems 
[13,14]. If we require the approximate parametric representation to give the correct (k — 1) 
universal ratios r 6 , r 8 , . . ., r 2 k+2, we obtain 

Wi = E c nm 6 m (h 3 + 7)" ^' , (4.6) 

m=0 l Zm + 

where 

c ™ = u—^v II ( 2 0™ - 7 + * - 1), (4-7) 

fc=i 

and we have set r 2 = r 4 = 1. Moreover, by requiring that F(z) = z + |-2 3 + we obtain 
the relation 

p 2 = 6(/i 3 + 7)- (4.8) 

In the exact parametric representation, the coefficient h 3 can be chosen arbitrarily. This is no 
longer true when we use our truncated function h(9), and the related approximate function 
/ipprox( x 5 ^3) depends on /i 3 . We must thus fix a particular value for this parameter. In order 
to optimize this choice, we employ a variational procedure [14], requiring the approximate 
function f^pl rox (x, h 3 ) to have the smallest possible dependence on h 3 . This is achieved by 
setting /i 3 = hz,ki where h^k is a solution of the global stationarity condition 
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TABLE I. Expansion coefficients for the scaling equation of state obtained by the k = 1,2 
approximations of the parametric function h{9). See text for definitions. The number marked with 
an asterisk is an input, not a prediction. 





k = 1 


k = 2 


Ising (Ref. [15]) 


/l° 


1.100(2) 


1.06(2) 


1.0527(7) 


f°2 


0.083(1) 


0.06(2) 


0.0446(4) 


fi 


-0.012(2) 


-0.006(3) 


-0.0059(2) 


ycoex 


0.87(1) 


0.93(8) 


0.9357(11) 


Jo 


0.550(3) 


0.48(3) 


0.6024(15) 


r% 


1.46(6) 


*0.90(15) 


2.056(5) 


T8 


1.0(3) 


1.5(3) 


2.3(1) 


V3 


6.9(2) 


6.3(2) 


6.050(13) 


V4 


17(1) 


18.7(5) 


16.17(10) 


jpoo 
^0 


0.0235(6) 


0.018(2) 


0.03382(15) 



dfW rox (x,h 3 ) 



= (4.9) 

^3=^3, k 



for all x. The existence of such a value of h 3t k for each A; is a nontrivial mathematical 
result which was proved in Ref. [14]. This procedure represents a systematic approximation 
scheme, which is only limited by the number of known terms in the small-magnetization 
expansion of the Helmholtz free energy. Note, for k — 1, the so-called linear model, Eq. 
(4.9) gives 

^« 

which was considered as the optimal value of h% for the Ising equation of state [37]. 



B. Results 



We apply the method outlined in the preceding section, using the Monte Carlo estimates 
[17] v = 0.683(3), r\ = 0.035(2), r 6 = 0.90(15) as input parameters. We obtain two different 
approximations corresponding to k — 1 and k = 2. Using the central values of the input 
parameters, we have 

h(6) ik=l ^ = e(i-e 2 /el), e 2 = 1.61451, (4.11) 

that provides the optimal linear model, and 

h(ef k= ^ = 9 (1 - e 2 /el) (1 + Cl e 2 ) , (4.12) 

9l = 1.32141, ci = 0.106612. (4.13) 

The relatively small value of C\ supports the effectiveness of the approximation scheme. 
In Table I we report results concerning the behavior of the scaling function f(x), F(z), 
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FIG. 1. The scaling function f(x). 

and $(u) for H = and on the critical isotherm, cf. Eqs. (2.3), (2.4), (2.5), (2.8), (2.9), 
(2.10). The errors reported there are only related to the uncertainty on the corresponding 
input parameters. We consider the k = 2 results as our best estimates. Of course, the 
corresponding errors do not take into account the systematic error due to the approximation 
scheme. Nevertheless, on the basis of the preceding applications [14-16] to the three- and 
two-dimensional Ising universality class, we believe that the k = 2 approximation already 
provides a reliable estimate, and that one may take the difference with the k — 1 result as 
indicative estimate of (or bound on) the systematic error. For comparison, the last column 
reports the corresponding estimates for the Ising universality class, taken from Ref. [15]. 
In Fig. 1 we show the scaling function f(x), as obtained from the k = 1,2 approximations 
of h(9), using the central values of the input parameters. The difference between the two 
curves is rather small in the region — 1 < x ^ 1. For larger values of x some differences 
are observed: they are essentially due to the small difference (approximately 10%) in the 
corresponding values of f£°. 

V. UNIVERSAL AMPLITUDE RATIOS 

Universal amplitude ratios characterize the critical behavior of thermodynamic quantities 
that do not depend on the normalizations of the external (magnetic) field, of the order 
parameter (magnetization), and of the temperature. From the scaling function f(x) one 
may derive many universal amplitude ratios involving quantities defined at zero momentum 
(i.e. integrated in the volume), such as the specific heat, the magnetic susceptibility, etc.... 
For example, one can obtain estimates for the specific-heat amplitude ratio A + /A~, the 
sueceptibility amplitude ratio C + /C~, etc.... Then, using the results for G 4 and the relation 
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TABLE II. Universal amplitude ratios obtained from the k = 1, 2 approximations of the para- 
metric function h(6). Amplitude definitions are reported in App. A. For comparison, the last 
column reports the corresponding estimates for the Ising universality class, taken from Ref. [15]. 



k = 1 k = 2 Ising 

A+/A- 1.33(8) 1.58(26) 0.532(3) 

C+/C- 4.5(2) 5.5(5) 4.76(2) 

R+ = aA + C + /B 2 0.1008(5) 0.079(6) 0.0567(3) 

#4 = -C+B 2 /(C+) 3 9.9(2) 12.2(8) 7.81(2) 

R x = C + B s - l /(B c ) 5 1.82(1) 2.1(1) 1.660(4) 

G\ /3 R+ = (-aA+C+/(C + ) 2 ) 1 / 3 0.999(5) 0.989(5) 0.443(2) 



ct 

° 4 " ~(C+) 2 (/+) 3 ' ^'^ 

where C4", C + , and / + are respectively the amplitudes of the zero-momentum connected 
four-point correlation, the susceptibility, and the second-moment correlation length (see 
App. A for notations), we can also determine universal ratios involving the correlation- 
length amplitude /+. 

In Table II we report estimates of several amplitude ratios, as derived by using the ap- 
proximate polynomial representations of the equation of state for k = 1,2. Again, the errors 
reported there are only related to the uncertainty on the corresponding input parameters. 
Note that, in the most important case of A + /A~, the k = 2 estimate includes the k = 1 
result within its error. As discussed in Sec. IV B, we consider the k — 2 result as our best 
estimate for each quantity. As an indicative estimate of the total uncertainty, we take the 
maximum between the error induced by the input parameters and the difference with the 
k = 1 result. This would lead to the final estimates A + / A~ = 1.6(3), C + /C~ = 5.5(1.0), 
R+ = 0.08(3), Rt = 12(2), R x = 2.1(3) and Gl /3 R£ = 0.99(1). 

Universal amplitude ratios involving the correlation-length amplitude / + , such as 

R+ = M+) 1 / 3 / + , Qc = B 2 {f+f/C + , (5.2) 
can be obtained by using the Monte Carlo estimate of Ref. [17], G 4 = 43.3(2): 

Rt = \^r-) = °- 282 ( 3 )' ( 5 - 3 ) 

Q c = S- = 0.28(5). (5.4) 

The result for is in substantial agreement with the rather precise Monte Carlo es- 
timate [17] R£ = 0.2885(15), providing support to our approximation of the equation of 
state. We have also obtained another independent estimate of the ratio R£ by using the FT 
fixed-dimension approach in terms of zero-momentum renormalized coulings. As discussed 
in detail in App. C we have extended the five-loop calculation of Ref. [12] to six loops. The 
analysis of the perturbative expansions gives the estimate R£ = 0.290(10), which is in good 
agreement with the results from the equation of state and from Monte Carlo simulations. 
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The comparison with experiments on uniaxial antiferromagnets, see, e.g., Ref. [3] for a 
review, is essentially restricted to the specific-heat ratio A + /A~ . Our best estimate A + / A~ = 
1.6(3) is in good agreement with the experimental result A + /A~ = 1.6(3) reported in Ref. [7], 
and A + /A~ = 1.55(15) reported in Ref. [3]. Earlier FT studies in the framework of e and 
fixed-dimension expansions [9-11] provided rather imprecise results. Actually, they favored 
a negative value, in substantial disagreement with experiments. 

Other experimental results concern the staggered magnetic susceptibility and the correla- 
tion length, see, e.g., Refs. [38,3], which are determined from neutron-scattering experiments. 
As noticed in Ref. [39], the two-point correlation function measured in these experiments 
does not coincide with 

G(x) = (p s p x s x ) c = (p s p x s x ) - ^ piSi) 2 , (5.5) 

i 

whose zero-momentum component is given by the susceptibility \ = dM/dH, but rather 
with 

1 -— 2 

G(x) = (p S pxS x } - y(22pi S i) ■ ( 5 - 6 ) 

i 

This difference affects essentially the low-temperature critical behavior. Indeed, setting 
X = C ± \t\~' y for the critical behavior of the corresponding susceptibility in the high- and 
low-temperature phase, C + = C + , but C~ ^ C~ , and therefore C + /C~ ^ C + /C~. A 
leading-order calculation [39] within the ^fe expansion gives 

(5.7) 

Trusting the leading-order result (5.7) and using C + /C~ = 5.5(1.0), see Table II, we obtain 
approximately C + /C~ w 3, which is consistent with the experimental result C + /C~ = 2.8(2) 
of Ref. [38]. 



c- 



APPENDIX A: NOTATIONS 

In this appendix we define the thermodynamic quantities considered in this paper and 
their critical behavior. We consider: the specific heat 

Ch = \ (M 5 ) , (Al) 

(E = —JJ2<ij> Pi Pj SiSj), whose critical behavior is 

C H = D + A ± \t\~ a , (A2) 

where D is a background constant that is the leading contribution for t — > since a < 0; 
the spontaneous magnetization near the coexistence curve 

M = ±-(EPiSi) = B\t\-e; (A3) 
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the magnetic susceptibility x an d the second-moment correlation length £, 

x = C ± \t\-\ t = f ± \t\- v , (A4) 
defined from the connected two-point function averaged over random dilution 



G{x) = (pos p x s x ) c ; (A5) 

the n-point susceptibilities Xn, defined as the zero-momentum connected n-point correlations 
averaged over random dilution, whose asymptotic critical behavior is written as 



Xn 



= C^r 7 ~ ( "~ 2)/M (A6) 



We also consider amplitudes defined in terms of the critical behavior along the critical 
isotherm t — 0, such as 

M = B C H 1/S , X =^\H\- J/I3S . (A7) 

o 



APPENDIX B: CORRESPONDENCE BETWEEN RIM AND FT 

CORRELATIONS 

In this appendix we determine the relationships between the RIM correlation functions of 
the spins and those of the field (p(x) that can be computed in the FT approach discussed in 
Sec. Ill B. In the RIM, beside the n-point susceptibilities Xn, defined as the zero-momentum 
n-point connected correlation functions averaged over random dilution, one may also con- 
sider dilution averages of products of sample averages. 

Setting 

**^(I>0*) , ( B1 ) 

we define the dilution-averaged correlations 

Pk!k 2 ...k n = cr kl a k2 ...a kn , (B2) 

and the generalized susceptibilities 
1 

Xk = yPki 

Xk x k 2 = y (Pk!k 2 - Pk 1 Pk 2 ) , (B3) 
Xk x k 2 k z = y (Pfcifefes - PkxkiPkz ~ Pkik 3 Pk 2 ~ Pk 2 k 3 Pki + ^PkxPk 2 Pk z ) , 

etc.... The fc-point susceptibilities Xk have already been introduced in App. A. Analogously, 
beside the universal quantities G k defined in Eqs. (3.4) we consider the quantities 
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G22 = lim lim 
G42 = lim lim 

t^0+ V^oo 



G 



222 



lim lim 

t^0+ V^oo 



X22 
X42 

'e d xi 
X222 



+4 



+ 6 



X4X22 

e d xi 

X22 



where V is the volume. In the FT approach one starts from the Hamiltonian 



d 3 x 



(B4) 



(B5) 



where ip(x) is a spatially uncorrelated random field with Gaussian distribution. In the limit 
of small dilution this FT model should have the same critical behavior of the RIM [26]. 
Setting 



ft 



d 3 x 4>(x) 



(B6) 



we define kn and Xkl,...,k n by using Eqs. (B2) and (B3), with a k replaced by cr^ T . Here the 
overline indicates the average over the random field ij)(x). The universal constants G^, ...,&„ 
are given by the same expressions used in the spin model, Eqs. (3.4) and (B4). 

If Z^(H) is the partition function (the generator of the connected correlation functions 
in the FT language) for given H and disorder configuration if), we define 



so that 



FT 

Pk x ...k 



A r (H l ,...,H r ) = lnZ^tfO ...InZ^Hr), 



-A r (Hi, ...H r 



#1=0,... #,.=0 



To compute these quantities we use the standard replica trick, i.e. we rewrite 

Z^{H) N -\ 



\nZj,(H) = lim 



N 



Introducing rN replicas, we write 



A r (Hx, ...H r ) 



1 



[dip] [dipij] . . . [d<f r;i \ exp 



N 



k=l i=l 



where 



Ti 



rip 



(B7) 



(B8) 



(B9) 



(BIO) 



(Bll) 



ki 



k = 1, ...r, and % — 1, ..N. In Eq. (BIO) we retained only the term depending on all arguments 
since we only use this expression to compute pj^ kr for ki > for all i. Integrating out the 
disorder field if), we obtain the generator of the connected correlation functions 
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A r (H 1 ,...H r ) oc J Y[[d(p k ,i\exp 



N 



-H rifi 4 + E E 

k=l i=l 



(B12) 



where 



7-V = / ^Ue [(d^) 2 + r^] + i«o (E + E ¥>« } • (B13) 

The Hamiltonian (B13) is equivalent to the Hamiltonian (3.9) with iV replaced by rN. In 
the limit N — > 0, this is of course irrelevant. Thus, if we define the n-point connected 
correlation function for the theory with Hamiltonian (3.9), 



*£..,a n = (<Pa 1 ---<Pa n ) c , 



we obtain 



N 

Y — ij m _L V X {n) 

ai,...,a n =l 
1 " 



27V 



lim 



1 



ai...a n bi...6 m ' 



E E 

oi,...,a„=l 6i,...,6 m =iV+l 

AT 2AT 3JV 

Xnmp = Jhn ^3 EE E v-v ai...a n bi...fe m ci...c p - 

oi,...,a„=l 6i,...,6 m =JV+l ci,...,c p =2JV+l 



(n+m+p) 



(B14) 

(B15) 
(B16) 
(B17) 



We wish finally to relate the constants G^...^ with the universal constants that parametrize 
the Helmoltz free energy of the FT model, cf. Eq. (3.10). Using 



x (n) 

a\...a n 



d n ~ l M n 



0H a2 ...dHj 



(B18) 



we obtain 



* { 2 = X2S a 



^abcd — S A 2 



94i- {Sabred + sym) + globed 



xSLf = -? d xl[(g e i - lO^^^cA/ + sym) + 

(962 - 20# 4 is> 42 ) -^((UcA/ + sym) + (g 63 - I0gl 2 )5 abcdef 



1 

15 



(B19) 



where S ai ,„ an = 1 if all indices are equal and zero otherwise, and "sym" indicates the appro- 
priate permutations. 

Then, using Eqs. (3.4) and (B4), we find the following relations 
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G 4 = 342 = v*, (B20) 

G22 = ^4i = \u\ (B21) 

G 6 = (? 63 , (B22) 

G42 = ^r<?62, (B23) 

G222 = ^61, (B24) 
G s = 9 85 - (B25) 

Let us summarize the available numerical estimates for the above-reported quantities. The 
Monte Carlo simulations reported in Ref. [17] provided the results: G 4 = 43.3(2), G22 = 
-6.1(1), r 6 = 0.90(15), C 42 = G 42 /(G 4 G 22 ) = 0.12(5), and C 222 = G222/G22 = 0.45(15). For 
the six-point couplings, we obtain correspondigly Gq = 1.7(3) x 10 3 , G42 = —32(13), G222 = 
17(6). The available FT estimates are more imprecise. For the four-point couplings, Ref. [19] 
applied different resummation methods to the six-loop /5-functions, obtaining estimates that 
can be summarized by G 4 = 38.0(1.5), G22 = —4.5(6). These estimates differ significantly 
from the Monte Carlo ones, a discrepancy that is probably due to the non-Borel summability 
of the perturbative series [40]. Estimates of Gq can be obtained from the results for r 6 
reported in Sec. Ill B. By using the Monte Carlo estimate of G4, we obtain Gq = 2100tgoo. 
Results less dependent on the location of the fixed point can be obtained by multiplying 
the estimate of r 6 at the FT (resp. Monte Carlo) fixed point for the corresponding FT 
(resp. Monte Carlo) estimate of G 4 . At the FT fixed point G 4 « 38 and r 6 = 1.1(1) so that 
G 6 = 1.6(2) x 10 3 , while at the Monte Carlo fixed point G 4 = 43.3(2) and r 6 = 0.6(3) that 
implies Gq = 1-1(5) x 10 3 . Taking as final estimate that derived by using the FT estimate 
of G 4 , we have G e = 1600i5oo, where the error is such to include also the estimate at the 
Monte Carlo fixed point. We also used field theory to determine the replica-replica six-point 
couplings. We applied the Pade-Borel method to the four-loop perturbative expansions 
[28,29] of C 42 = ^42/(^22^4) and of C222 = G222IGW. The results are only indicative since 
the estimates change significantly with the order and with the Pade approximant. We obtain 
C 42 = 0.06(3), 0.01(5) and C 22 2 = 0.05(2), -0.015(25) by using the FT and the Monte Carlo 
estimate of the fixed point respectively. Comparing with the Monte Carlo results reported 
above we observe that field theory (at four loops) provides only the order of magnitude, 
but is unable to be quantitatively predictive. Finally, we consider r 8 and Gg. As discussed 
in Sec. Ill B (see also Ref. [35]), the perturbative three- loop expansion of r 8 gives only an 
upper bound. A much more precise estimate of r 8 has been derived in Sec. IVB, r 8 = 1.5(5), 
which implies G 8 = 1-2(4) x 10 5 . 



APPENDIX C: FIELD-THEORETICAL DETERMINATION OF i?+ 

In this appendix we estimate the universal amplitude ratio = f + (aA + ) 1 ^ 3 by per- 
forming a six-loop expansion in the framework of a fixed- dimension FT approach based on 
the Hamiltonian (3.9). In this scheme the theory is renormalized by introducing a set of zero- 
momentum conditions for the one-particle irreducible two-point and four-point correlation 
functions: 
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r { 2{p) = s ab zz 1 U 2 + p 2 + o(p 4 ) 



ritL(o) 



zr 2 



m (uS abcd + vC abcd ) 



(CI) 
(C2) 



where S abcd = | (S ab 5 cd + 5 ac 5 bd + S ad 5 bc ) and C abcd = S ab 5 ac 5 ad . Eqs. (CI) and (C2) relate 
the mass scale m and the zero-momentum quartic couplings u and v to the corresponding 
Hamiltonian parameters r, u , and v , 



uq = muZ u Z~ 



v = mvZ v Z' 



(C3) 



In addition one defines the function Z t through the relation 

rii 2) (o) = W, 



(C4) 



where r^ 1 ' 2 ) is the one-particle irreducible two-point function with an insertion of |0 2 . The 
pertubative expansions of Z v (u,v), Z u (u,v), Z v (u,v), Z t (u,v) have been computed to six 
loops [27,19]. The specific heat is given by the zero-momentum energy-energy correlation 
function averaged over the random dilution. In the FT approach it corresponds to 



Cm = lim 



1 



n^o N 



d d x (^(OV 2 (rr)> ( 



(C5) 



i.e. to the zero-momentum value of the two-point correlation function of the operator (p 2 = 
Yja^a- We computed Ch to six loops, extending the five-loop computation of Ref. [12]. 
The calculation requires the evaluation of a few hundred Feynman diagrams. We handled 
it with a symbolic manipulation program, which generates the diagrams and computes the 
symmetry and group factors of each of them. We used the numerical results compiled in 
Ref. [41] for the integrals associated with each diagram. 

In order to compute the universal ratio R^, we follow Ref. [11]. We consider different 
expansions that converge to as t — > + : 



R 



(i) 



lim 



ia-c H ) i/d 



dt 



(C6) 
(C7) 



where C' H = dCn/dt (at fixed uq and vq) and t is the reduced temperature. The derivatives 
with respect to the reduced temperature t can be done by using Eq. (C4), which can be 
rewritten as 



dt 



zr 1 - 



(C8) 



u ,v 



This allows us to compute the derivative with respect to t of generic functions written in 
terms of Uq,vq, and x- For example, setting 



C H = x 1/2 Ch(u 1 , vi), ui 



X 1/2 u , 



Vi 



X 1/2 vo, 



(C9) 
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we obtain 



Cff(t) = -\x V2 z; 1 [1 + Ul ^- + «i^J ChK,^i). (CIO) 

Using the relations £ = 1/m, x — -Z^£ 2 , an d Eq. (C3), one obtains expressions that can be 
expanded in powers of the renormalized quartic couplings u and v. For example, 

Rf\u,v) = (l + +vi-^j Ch^v,), (Cll) 

where u± = Z u Z~ l u, v\ = Z v Z^ p 1 v. We write 

<W) = Ei n) *M (C12) 

where u and v are the rescaled couplings u = uj (Qir) and v = 3t>/ (167r). The coefficients cj™' 1 
are reported in Table III for < i+j < 5. Estimates of R£ are obtained by resumming these 
series, and evaluating them at u* and v*. We employed several resummation methods, see, 
e.g., Refs. [19,17]. Without entering into the details, we report the results: = 0.287(2) 

and = 0.286(3) obtained by using the FT estimates of the fixed point, u* = 37 and 

v* = -13, and Rf = 0.296(2) and Rf = 0.295(2) obtained by using the Monte Carlo 
estimates [17] u* = -18.6(3) and v* = 43.3(2). Our final estimate is R^ = 0.290(10) that 
includes all the above-reported results. Our FT estimate agrees with the Monte Carlo result 
of Ref. [17], R£ = 0.2885(15), and with the FT estimates of Refs. [11,12], R+ = 0.286(4) 
(four loops) and R^ 0.2887 (five loops). 
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TABLE III. Coefficients cf/ and &^ of the expansion of R\ and R\ respectively. 



i,j 


(i) 
r. ■ 


(2) 


0,0 


0.2150635 


0.2150635 


0,1 


0.0358439 


0.0358439 


1,0 


0.0268829 


0.0268829 


0,2 


0.0000532 


-0.0019913 


1,1 


-0.0043608 


-0.0089610 


2,0 


-0.0016353 


-0.0033604 


0,3 


0.0025329 


0.0055607 


1,2 


0.0079130 


0.0183873 


2,1 


0.0079586 


0.0187523 


3,0 


0.0019897 


0.0046881 


0,4 


-0.0021046 


-0.0052579 


1,3 


-0.0097350 


-0.0242605 


2,2 


-0.0167427 


-0.0414335 


3,1 


-0.0110331 


-0.0273141 


4,0 


-0.0020687 


-0.0051214 


0,5 


0.0024689 


0.0072264 


1,4 


0.0140199 


0.0411159 


2,3 


0.0316913 


0.0929392 


3,2 


0.0345501 


0.1009750 


4,1 


0.0165332 


0.0482474 


5,0 


0.0024800 


0.0072371 
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